function [params_out,end_moments,best_distance] = sim_annealing_parallel(func,targets,x0,ub,lb,max_loops)

    start_t=now;
    num_targets=size(targets,1);
    num_top_params=4;
    num_params=size(x0,1);
    start_x_rep=nan(num_params,num_top_params^2);
    start_x=x0;
   
    distance=nan(num_top_params^2 + num_top_params,1);
    
    start_targets=nan(num_targets,num_top_params);
    attempt_targets=nan(num_targets,num_top_params^2 + num_top_params);
    target_mat=repmat(targets,1,num_top_params^2);
    parfor y=1:num_top_params
        disp(y)
        [distance(y),start_targets(:,y)]=func(x0(:,y),target_mat(:,y));
    end
    %tmpdistance=nan(num_top_params^2 + num_top_params,1);
    %x_step = nan(num_params,num_top_parms^2);
    %x_step_pos = nan(num_params,num_top_parms^2);
    %x_step_neg = nan(num_params,num_top_parms^2);
    for y=1:num_top_params        
        start_x_rep(:,num_top_params*(y-1)+1:y*num_top_params)= repmat(x0(:,y),1,num_top_params);
    end
    
    for j=1:max_loops
        loop_start=now;
        parfor z=1:num_top_params^2
            disp(z)
            rng(sum(100*clock*z))
            x_step = ((ub-lb).*rand(length(lb),1)-(0.5*(ub-lb)))*(sqrt((max_loops-j)/max_loops));
            on_lb = start_x_rep(:,z)==lb;
            x_step_pos = ((ub-lb).*rand(length(lb),1))*(sqrt((max_loops-j)/max_loops));
            x_step(on_lb)=x_step_pos(on_lb);
            on_ub = start_x_rep(:,z)==ub;
            x_step_neg = ((lb-ub).*rand(length(lb),1))*(sqrt((max_loops-j)/max_loops));
            x_step(on_ub)=x_step_neg(on_ub);
            disp(['z=' num2str(z) ',xstep is:'])
            disp(x_step)
            x_guess(:,z) = min(max(start_x_rep(:,z) + x_step,lb),ub);
            [distance(z+num_top_params),attempt_targets(:,z+num_top_params)] = func(x_guess(:,z),target_mat(:,z));
        end
        loop_end=now;
%         if (loop_end-loop_start)*1440>5
%             mail_text=['theta val=' num2str(x0(1,end)) ', time for this loop=' num2str((loop_end-loop_start)*1440)];
%             system(['mail -s "Each sim annealing parallel loop is going really slowly" inath33@gmail.com<<<' mail_text])
%         end
        [~,best_choices]=sort(distance,1,'ascend');
        attempt_targets(:,1:num_top_params)=start_targets;
        full_x_guess=[start_x,x_guess];
        for zz=1:num_top_params
            start_x_rep(:,num_top_params*(zz-1)+1:zz*num_top_params)=repmat(full_x_guess(:,best_choices(zz)),1,num_top_params);
            start_x(:,zz)=full_x_guess(:,best_choices(zz));
            tmpdistance(zz)=distance(best_choices(zz));
            start_targets(:,zz)=attempt_targets(:,best_choices(zz));
        end
            distance(1:num_top_params)=tmpdistance;
            disp('')
            disp(['Loop num=' num2str(j) ])
            disp('Best Distances:')
            disp(distance(1:num_top_params)')
            disp('Best Start params:')
            disp(start_x(:,1:num_top_params))
            disp('Best Start targets:')
            disp(start_targets(:,1:num_top_params))
 
    end
        params_out=start_x_rep;
        end_moments = start_targets(:,1);
        best_distance = distance(1);
        disp(['total minutes = ' num2str((now -start_t)*(86400/60))])
    end
    